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Abstract 

Modified gravity theories with the Gauss-Bonnet term G = B? — AR^^^R^u + R^'^p^ R^^pcr have 
recently gained a lot of attention as a possible explanation of dark energy. We perform a thorough 
phase space analysis on the so-called f{G) models, where f{G) is some general function of the 
Gauss-Bonnet term, and derive conditions for the cosmological viability of /(G) dark energy mod- 
els. Following the f{R) case, we show that these conditions can be nicely presented as geometrical 
constraints on the derivatives of /(G). We find that for general /(G) models there are two kinds of 
stable accelerated solutions, a de Sitter solution and a phantom-like solution. They co-exist with 
each other and which solution the universe evolves to depends on the initial conditions. Finally, 
several toy models of /(G) dark energy are explored. Cosmologically viable trajectories that mimic 
O the ACDM model in the radiation and matter dominated periods, but have distinctive signatures 

^ at late times, are obtained. 
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I. INTRODUCTION 



The discovery that our universe recently started accelerating PQ [21 El H] has baffled 
particle physicists and cosmologists for a decade. Many suggestions about the origin of the 
cosmic acceleration, or the nature of dark energy, have been proposed (for a review, see [S]). 
While the cosmological constant (ACDM) introduced by Einstein is the simplest proposal, 
albeit requiring fine-tuning to avoid coincidence problems other exotic negative-pressure 
fluids, often described in terms of scalar fields [7j, have been proposed to address these 
issues. However, it is probably fair to say that no compelling and well-motivated solutions 
have yet been developed. As a result, attention has naturally turned to the possibility of 
dark energy originating from modifications of gravity itself [H [9l [TOl [TTl [121 [13] . 

Among the recent modifications of general relativity, perhaps the simplest approach is 
to replace the Einstein-Hilbert Lagrangian density with some general function of the Ricci 
scaler v^/(i?) [ni[ini[ia[I51[ISl[ni[ISl[I2]. For these kind of models, generally it is easy to 
obtain a late time accelerated epoch. In [2ni[2T], conditions for the cosmological viability of 
general f{R) models have been established. But as with the popular model f{R) = R+fi/R, 
they are also generally found to be inconsistent with solar system constraints [221 ES] . 

A well motivated curvature invariant beyond the Ricci scalar is the the Gauss-Bonnet 
term 

G = R'- AR^'R^, + R^'^'^R^up. , (1) 

which is inspired by string theory [211 [2S] and has gained special interest in cosmology 
[2S1 [23 [2H1 [23 [201 [SI] • it is also well known that the Gauss-Bonnet combination can reduce 
the number of spin-2 ghosts, which otherwise severely haunt the perturbation theory [32]- 
The 4 dimensional Gauss-Bonnet term is a topological invariant, thus has no dynamical 
effects if added into the Lagrangian linearly. To introduce extra dynamics, one may couple 
the Gauss-Bonnet term to a scalar field, as naturally appears in low-energy string effective 
actions [24] . For the exponential coupling with a scalar field potential, this model can 
produce a matter dominated period followed by an accelerated period [33l|3l]. However, it 
has been shown that it tends to conflict with solar system constraints [35] . 

Alternatively, we may consider the Lagrangian density as some general function of the 
Ricci scalar and the Gauss-Bonnet term y/—gF{R, G) [121 [31], and as a natural simplifica- 
tion, F{R,G) = R/2 + f{G) has been investigated [301 [Ml [371 [Ml [39]. Uddin et al. [10] 
studied the existence and stability of power-law scaling solutions in f{G) models, finding 
that scaling solutions exist in the model f{G) = ±2V aG with a an arbitrary constant. 
In [H], the authors derived the stability conditions for a de Sitter solution and a stan- 
dard matter/radiation solution in general f{G) models; the authors also constructed several 
cosmo logically viable f{G) models, but found it difflcult to simulate the models to high 
redshifts. 

In this paper, we systematically survey the cosmological viability of general f{G) models. 
The layout of this paper is as follows. In Section [Tl| we first derive field equations of general 
f{G) models in the fiat Robertson- Walker background and recast them nicely into a 5 
dimensional autonomous system. Analogous to the case of f{R) models [20], we construct 
the curve m(r), where m = Gfccl fc and r = —Gfc/f with fc = d//dG, etc.. In Section 



III, the properties of the phase space and their cosmological implications are determined. 



Interestingly, instead of critical points, we obtain lines of critical points in the phase space. 
We then establish conditions for cosmo logically viable f{G) dark energy models in the m-r 
plane. Consequently, when deciding the cosmological viability of a particular f{G) model, 



2 



one's first resort may be to inspect its m(r) curve in the m-r plane. The particular results of 
can be easily recovered from the more general conditions obtained with our approach. 



In Section IV, we investigate several specific f{G) models and examine their cosmological 
viability by numerical simulations. Cosmologically viable trajectories of the universe are 
obtained. We summarise the results in Section |Vl 



II. DYNAMICAL FRAMEWORK FOR f{G) COSMOLOGY 
We consider the following action 

(2) 



S = J d^x^g Qi? + /(G) + Cr + 



where R is the Ricci scalar, f{G) is a general function of the Gauss-Bonnet term G and 
y/—gCr and \/—gCm are the radiation and matter Lagrangian densities respectively. We 
have chosen SuGh = 1, where Gb is a bare gravitational constant. Varying the action ^ 
with respect to g^u, we obtain the field equations for f{G) gravity 

Rfii/ '^^Rg^iv 

-|-8 [-R^p[y(j -|- Rpugafi Rpcrgup Rfiugap ~l~ Rfia^gup 

+ \R{9pug.p - 9p.9up)]^''V''fG + (G/g - f)9pu = T^u , (3) 

where T^j^ is the energy momentum tensor of radiation (energy density pr) and matter 
(energy density pm)- In the observationally favoured flat Robert son- Walker metric [1] 

ds^ = -dt^ + a(t)2dx2 , (4) 

we have 

R = 6{H + 2H^), G = 24H\H + H'^) , (5) 

and Eqs. ^ reduce to 

3H' = Gfc-f- 2AH'fG + Pr + Pm, (6) 

-2H = -SH^fc + IQHHfc + SH^fc + + Pm , (7) 

where H is the Hubble parameter and an overdot stands for a derivative with respect to t. 
Additionally, the densities pr and pm satisfy the usual continuity equations 

Pr + 4Hpr = 0, (8) 
Pm + SHpm = . (9) 

Eqs. ([6]|9| determine the dynamics of the f{G) gravity system ([2]) in a homogeneous and 
isotropic background. 

Drawing the analogy with 3H^ = p and —2H = p + p, we can naturally define gravita- 
tionally induced forms of dark energy density pde and pressure as 

PDE = GfG~f-24H'fG, (10) 

Pde = IQH^fc + IQHHfc + SH^/g -Gfc + f. (11) 
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In this way, the continuity equation of dark energy holds, 

pDE + ^H{pde + Pde) = , (12) 
and its equation of state parameter becomes 

Pde IGH^fc + 16M/g + SH^fg -Cfg + f 

wde = = -■ • (13) 

Pde GfG-f-24H^fG 

We also define the effective equation of state [20] 



and dimensionless energy densities 



2H , , 

^e// = -l-^, (14) 



f^x = ^, X = r,m,DE. (15) 

To discuss the dynamics of a general /(G) model, it is convenient to rewrite the equations 
of motion as a dynamical system jl2| HSj 01]. To achieve this we introduce the following 
dimensionless variables 

Gfc /,„x 

^1 = ^ , (16) 

^2 = -^, (17) 

X3 = -8iJ/G, (18) 

= ^r = 7^ , (19) 

Then Eqs. ([6||9| can be recast as a first order dynamical system 

dxi X3X5 



dA^ m 
dx2 _ a;3a;5 

dA^ m 

dA^ 
dx4 

diV 

dX5 XsX^ 



X3X5 - 2x1X5 + 2xi , (21) 
- 2x2X5 + 2x2 , (22) 



— X3 + 2x5 — X3X5 + 1 — 3xi — 3x2 + 2:4 , (23) 
-2x4 - 2x4X5 , (24) 
4x^ + 4x5 , (25) 



dA^ X\ra 

where A^ = ln(a/aj) (oj is the initial value of the scalar factor) and 



m = , (26) 

Jg 



Gfc xi 



f X2 

4 



(27) 



Note that Eq. ^ is recast as 



Xi — X2 — 2^3 — 2J4 



In terms of Xi, we can rewrite wde and Weff as 

-2x5- 



Wde 



^eff 



X4 



3(xi + X2+ X3) 

-^(2x5 + 1). 



(28) 

(29) 
(30) 



From the r.h.s. of Eq. (24), we can factor out X4, which suggests that X4 = is an 
invariant submanifold of the dynamical system, meaning the system can not go through 
the subspace X4 = and can only approach it asymptotically. We can also factor out x^ 
from the r.h.s. of Eq. (25). However, X5 = is not an invariant submanifold, because when 
X5 = (i.e., G = 0), the first term {—x^xl/xim = G/2AH^) on the r.h.s. of Eq. (25) is not 
necessarily zero. 

From Eq. ([27)) we can express G as a function of Xi/x2, and then substituting it into 



Eq. (26), m can be expressed in terms of Xi/x2- Therefore, given a form of f{G), we obtain 
a function m(r = X1/X2) and the dynamical system (21 25| becomes autonomous. For 
example, the model f{G) = a{G'^ — f3Y corresponds to the straight line m(r) = (1 — q)r/q + 
p — 1 in the m-r plane. A similar function of m(r) exists in the f{R) models and examining 
the cosmological viability of f{R) dark energy models according to the corresponding m(r) 
curve in the m-r plane has proven to be a fruitful approach [20] ■ On the other hand, if 
we know a form of m(r) that is cosmologically viable, one may use it to obtain f{G) from 
Eqs. (|26| and (|27|). 

From Eqs. (22) and ([25]), given a form of f{G), we may express G and H in terms of X2 



and X5, and substituting it into Eq. (21 ), we may then express xi in terms of X2 and X5, so xi 
is actually not an independent variable. However, beyond the power-law form f{G) = aG^, 
which allows Xi to be eliminated directly, even simple forms of f{G) yield a complicated 
dependence for xi on X2 and Xs, which is often impossible to obtain in closed form. Note 
that in [31], the authors recast the field equations into a 4 dimensional dynamical system 
(autonomous if the form of f{G) is given). However, phase space analysis, especially for 
general f{G) models, in that reformulation is not easy. Consequently we solve the original 5 



dimensional system (21 25), which is analytically simple and numerically more reliable and 
efficient. 



III. ANALYTICAL RESULTS OF GENERAL f{G) MODELS 



Writing the field equations in terms of the new dimensionless variables Xi makes it easy to 
look at the analytical properties of general f{G) models. The critical points (equilibria) are 
where the r.h.s. of the first order dynamical system (21 25) is zero. A new feature emerges 
here. Rather than having isolated critical points, we obtain continuous lines of critical points 
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(equilibrium manifolds), which we call critical lines: 



Li : {xi = 1 - X2, X2 = X2, X3 = 0, X4 = 0, X5 = 1} , 

fim = 0, ilr = 0, n£,£; = l, WdE = -1, Weff = -1, 
11 11, 

L2 ■ {xi = -X3, X2 = --X3, X3 = x3,Xi = 0,X5 = --,m = --} , 
bo II 

5 5 

^m = l--a;3, r^r = 0, ^DE=-X3-, WdE = , Weff=0, 
D D 

r r ^5 2X5 2(X5 - 1) 1 

L3 ■ {xi = -,X2 = -, X3 = — , X4 = 0, X5 = X5, m = --} 

X5 — 2 x^ — 2 X5 — 2 2 



ilm = 0, Qr = 0, fiz)£; = l, W DE = - 1^X5 - - , Weff = --X5 



2 1 _ 2 

11 3 1, 

L4 : {xi = -X3, X2 = "2^3, 2:3 = X3, X4 = 1 - -X3, X5 = -1, m = --} , 

3 3 11 

fim = 0, ^lr = l--X3, nDE=-X3, WdE=-, 1^eff = -- 

Li, L2 and L4 are straight lines in the phase space, while L3 is not. Since Li, L2 and L4 have 
different constant values of X5, 1, —1/2 and —1 respectively, they do not intersect. However, 
L3 intersects with Li at (—1, 2, 0, 0, 1), with L2 at (1/5, —2/5, 6/5, 0, —1/2) and with L4 at 
(1/3,-2/3,4/3,0,-1). 

Generally, if a nonlinear system has a critical line, the Jacobian matrix of the linearised 
system at a critical point on the line has a zero eigenvalue with an associated eigenvector 
tangent to the critical line at the chosen point. This kind of nonlinear system is a special 
subclass of the non-hyperbolic system (whose linearised system has one or more eigenvalues 
with zero real parts.). The emergence of equilibrium manifolds may be due to some symmetry 
in the nonlinear system or that the nonlinear system can be reduced to a lower dimensional 



one. For the case of f{G) = aG^, as we shall see in Section IV, when the system is reduced 
to a lower 4 dimensional one with just X2-X5, Li shrinks to a point but the other 3 lines 
remain. The stability of a particular critical point on the line can be determined by the non- 
zero eigenvalues, because near this critical point there is essentially no dynamics along the 
critical line (i.e., along the direction of the eigenvector associated with the zero eigenvalue), 
so the dynamics near this critical point may be viewed in a reduced phase space obtained by 
suppressing the zero eigenvalue direction. It is also interesting to note that the dynamical 



system (21 25) does not have any parameters in it, but, as we shall see below, the stability 
of a critical point changes along the critical line. Thus this system might be considered as 
the case of bifurcation without parameters 

Note that for L2, L3 and L4 to exist, m should be equal to —1/2, as we have noted in the 
definitions of the critical lines above. Also, for L2, L3 and L4, we have r = Xi/x2 = —1/2. 
Hence for these critical lines, 

M-^) = -^ (31) 

is satisfied, i.e., they must be located at the point (—1/2, —1/2) in the m-r plane. Thus for 
L2, L3 and L4 to exist, r = —1/2 should solve the equation m(r) = —1/2. From the values 
of floE, wde, etc., we can see that L2 and L4 correspond to solutions in which dark energy 
scales with matter and radiation respectively, thus they can hopefully correspond to the 
standard matter and radiation epochs. Therefore, for a standard matter or radiation epoch 
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to exist, the constraint (31) should be satisfied, i.e., the m(r) curve should pass through 
(—1/2,-1/2) in the m-r plane. 

Following the case in general f{R) models [20j, we can derive the following equation from 



Eqs. (26) and (27) 



_^ = r(m + r + l)^. (32) 

When r = is satisfied and (m + r + 1)G/HG does not diverge or m{r) = —r — 1 is 
satisfied and rG/HG does not diverge, we have dr/dN = 0. From the m-r plane's view, 
this means that, evolving along the curve m{r), typically the system can not go through 
any intersection points between the curve m(r) and the m axis (r = 0) or the particular 
straight line mc(r) = —r — 1. However, the curve m(r) can be "connected" at r = ±oo or 
m = ±oo. So for example, the system may evolve from the r < half plane to the r > 
half plane by "tunnelling" through the infinite point of the m(r) curve at r = ±oo. This 
happens because basically r is not an essential dynamical variable of the dynamical system 



(21 25). For example, r will go through the infinite point if X2 goes from 0~ to 0"^ with Xi 
remaining finite. 

Let us now discuss the properties and cosmological consequences of an arbitrary point 
on each of the critical lines in turn. Note that supposing the system starts from somewhere 
near a particular critical line and the evolution of the system is initially attracted to the 
critical line, then the particular critical point on the line the system will finally evolve to or 
pass nearby depends on the initial conditions. We denote an arbitrary point on, for example, 
Li as Pi. 

• Pi: (1 — 0:20, a;205 0, 0, 1), de Sitter point 

This is a de Sitter point for any point along the line Li, i.e., for any value of 0:20 (a 
given value of X2)- The eigenvalues of the linearised system are: 



where mi = m(ri = (1 — X2o)/3^2o)- Hence this is a stable spiral in the subspace of the 
last two eigenvalues when 

8 8 
< "^1 < N (^20 < 1) or — < mi < {x2o > 1) , (34) 

ZO[i — X20) ZO[i — X20) 



a stable node when 



< mi < — ^ (0:20 < 1) (35) 



or 



25 1-0:20 21-0:20 



1 8 

< "^1 < TTFTl N (^20 > 1) , (36) 



2(1-0:20) 25(1 -X20) 



and an unstable point otherwise. In the m-r plane, the stable spiral condition (34) 



corresponds to the area enclosed by the r axis and the curve mds*{r) = 8(r + l)/25r. 



whilst the stable node condition (35) corresponds to the area enclosed by the curve 



fnds*{f) = 8(r + l)/25r and the curve mds{r) = (r + l)/2r. Note that the point 
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(—1/2, —1/2) is on the curve mds{r) = (r + l)/2r (the boundary of the stable area), 
so the stable area of Li is connected to L2, and L4 in the m-r plane, as seen in 
Fig. [TJ One can check that these conditions are consistent with those established for 



a stable de Sitter point in [H]. For example, from the condition (34), we can derive 
the condition for the existence of a stable spiral in 



< HffGciG,) < ^ , (37) 

where Hi and Gi are the Hubble parameter and the Gauss-Bonnet term at the de 
Sitter point respectively and we have used 

GifaciGi) GifciGi) 4 

"^1 = r N ' 1 - X20 = xio = — —2 — , Gi = 2AHi . (38) 



• P2: (|a:3o, — 1^30, Xao, 0, — |), scaling with matter point 

At this point, dark energy mimics the evolution of matter, and the density of dark 
energy scales with that of matter {flm/^DE = (6 — 5a;3o)/5x3o), whereas there is no 
radiation. The eigenvalues of the linearised system are: 

0. -1. 3K(4).l). (30) 

where m'(— 1/2) = dm(r)/dr|r=-i/2- Along this line, the stability of the last two 
eigenvalues changes at the points X30 = or X30 = 96/71, which we may call bifurcation 
points: the points with X30 < or X30 > 96/71 are stable, while other points are 
unstable. In order for P2 to be a standard matter era {flm — 1), a^so needs to be 
near 0, but 0:30 = is a singularity of the last two eigenvalues. It is interesting to 
note that the situation here is quite similar to that of the standard matter point in 
general f{R) models [20], where for P5 in their notation to be a standard matter era 
needs to be near 0, but is a singularity of two eigenvalues. For our current 
point P2, if X30 0"*", the last two eigenvalues diverge as —00 and +00. However, 
if 0:30 is some small (non-infinitesimal) positive value, the two eigenvalues are large 
and finite with opposite signs. If the system flows towards these kind of points, unless 
the initial conditions are extremely fine-tuned, the system will not remain around it 
for a reasonable long time. Since an elongated matter dominated epoch is needed 
for cosmic structure formation, we generally consider this unacceptable. On 

the other hand, if X30 0^, then the last two eigenvalues are complex with negative 
real parts, which means the two eigenvalues are stable. However, the imaginary parts 
diverge, and hence the system oscillates rapidly in the subspace associated with the 
two eigenvalues when it approaches the standard matter point. The frequency of the 
oscillation can be reduced if 0:30 is slightly less than 0, while the amplitude of the 
oscillation will not be very large if the initial conditions are set near the critical point 
in the associated subspace. The system should finally leave the standard matter point 
and enter an accelerated point, so we shall set 3(m'(— 1/2) + 1) to be greater than 
in order to make the matter point unstable, i.e., 

m'{-^)>-l. (40) 
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f(G) = a -fa + p -fcf 



f(G) = a VcT exp(/3 G) 
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FIG. 1: The m(r) curves (red dashed hnes) for several /(G) models. Since these /(G) models 
involving functions such as \/G, G^^'^ or ln(/3G), G in /(G) should be understood as |G|. Note that 
the matter (radiation) point P2 (P4) and the phantom- like point P3 are located at (—1/2, —1/2). To 
have a standard matter dominated epoch, the curve should pass through the point (—1/2, —1/2), 
whilst not accessing the point from the forbidden directions (illustrated by the small light blue area 
around the point). The light green and light yellow areas are the potentially stable node and stable 
spiral regions of the de Sitter point Pi respectively. Evolving along the curve m(r), typically the 
system can not go through any intersection points between the m axis (r = 0) and the particular 
straight line mc{r) = —r — 1 (black solid line). However, the system can evolve through infinite 
points of the m(r) curve, as described in the main text. 



In the m-r plane, it means that there are unacceptable directions of the m(r) curve 
approaching the point (—1/2,-1/2), see Fig. [!} Note that the properties and their 
cosmological implications obtained above assume that the system is exactly at a critical 
point. Since the critical point here for the standard matter epoch is not supposed to 
be totally stable, the system will just pass nearby the critical point rather than fall 
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onto it, so the properties and their cosmological imphcations when the system passing 
by the critical point should be just approximate to what are inferred from the critical 
point analysis. 

(^f^) ^S^?' 0' ^50), dark energy dominated point 

This is a dark energy dominated point with no radiation and no matter. It has a 
fascinating wealth of possibilities, mimicking radiation when X50 = — 1 and matter 
when X50 = —1/2; it is a quintessence-like point (—1 < wde < —1/3) if < X50 < 1 
and a phantom-like point {wde < —1) if a^5o > 1- The eigenvalues of the linearised 
system are: 

0, -2x50-1, -X50-2, -2x50-2, -2(x5o-l)(m'(-^) + l). (41) 
This is a stable point when either of the following two conditions is satisfied 

- ^ < X50 < 1 and m'(-^) < -1 , (42) 

or 

X50 > 1 and m'(--) > -1. (43) 



Note that the second condition (43) is consistent with the condition (40), which means 
that if we require an unstable direction for the system to leave the standard matter 
point P2, then P3 can be a stable phantom-like point, co-existing with the stable de 
Sitter point Pi. Which point the universe finally evolves to depends on the initial 



conditions, as we shall see in Section IV, Note that although dark energy with w < —1 



is problematic theoretically |17], it is still consistent with current data |^ H8| H9]. 
P4:- {\x30, "^2:30, X30, 1 — fa^so, —1), scaling with radiation point 

At this point, dark energy mimics the evolution of radiation, and the density of dark 
energy scales with that of radiation {flr/^DE = (4 — 3x3o)/3x3o), whereas there is no 
matter. The eigenvalues of the linearised system are: 



, 1, , 1 1 /64 - 47x30 . s 

0. 1. 4(mV2) + l), -2±2V^^- 

The stability structure in the subspace of the last two eigenvalues is similar to that of 
P2: the points with X30 < or X30 > 64/47 are stable in the subspace (still unstable 
in the total phase space); in order for P4 to be a standard radiation point {Qr — 
1), X30 needs to be near 0, but at the same time, X30 = is a singularity of the 
two eigenvalues; if X30 is slightly greater than 0, P2 is extremely unstable; if X30 is 
slightly less than 0, P2 is stable in the subspace but the system oscillates significantly 
when flowing to it. Note that the condition fee > from [41] for the existence of 
a standard radiation/matter point can also be recovered. For example, to avoid the 
violent instability of the standard radiation point, we require X30 be slightly less than 0, 
so —f/SH"^ = X20 = —X30/2 > 0, which means / < 0. Additionally, for the standard 
radiation point, we have r = —Gfc/f = —1/2 < and m(— 1/2) = Gfcal fa = 
— 1/2 < 0, so fcG > holds at the standard radiation point. 



10 



Starting from the radiation dominated era, a cosmologically viable trajectory of the uni- 
verse in the phase space would start somewhere near the standard radiation point P4 (with 
X30 shghtly less than 0), then slowly pass nearby the standard matter point P2 (with X30 
slightly less than 0) and finally land at the stable de Sitter point Pi or the stable phantom- 
like point P3. In the m-r plane, with the trajectory fixed along the curve m(r) and subject 
to the constraint from Eq. (32), the universe would start near the point (—1/2,-1/2) where 
the standard radiation point P4 and the standard matter point P2 are located, slowly moving 
away from it, and finally fall onto some point in the stable area of Pi or come back to evolve 
to P3 that also resides at the point (—1/2, —1/2). Note that it might also be possible that 
the universe has an unstable accelelated epoch after the matter dominated period, whose 
cosmological viability is not the concern of this paper and is left for future work. Indeed, 
for phase spaces higher than 2 dimensions, more complicated dynamical phenomena such as 
chaos may appear, but again this is beyond the scope of the current paper. 



IV. SPECIFIC MODELS AND NUMERICAL RESULTS 

Having described possible cosmic trajectories for general f{G) dark energy models, we 
now turn our attention to the cosmological viability of a few specific toy models of f{G) 
whose m(r) curves can be analytically obtained. Note that if the analytical form of m(r) 
can not be obtained, one may get a numerical expression of it for a certain range of r. On 
the other hand, given a form of m(r), it may not be easy to obtain an analytical form for 
f{G), but again we can solve it numerically. 

The Gauss-Bonnet term G evolves from negative values in the radiation and matter 
dominated epochs to positive values in the acceleration epoch, thus the lagrangian is not 
well defined if there is any function of the Gauss-Bonnet term such as ^/G or \n(f3G). In 
such cases, G in the f{G) should be understood as |G| in this section as well as in Fig.Q We 
will start the system deep in the radiation dominated epoch, which requires we initially set 
r = Xi/x2 — —1/2, X4 ~ 1, ~ —1 and X1 + X2 + X3 ~ 0. Since the matter dominated point 
is located at (—1/2, —1/2), setting r = Xi/x2 — —1/2 is also a requirement for the universe 
to have an elongated matter dominated period. We find that the dynamical systems of f{G) 
models are quite often stiff. The simulation figures in this section are produced using the 
Matlab stiff ordinary differential equation solver odelSs. 

We first consider the simplest case f{G) = aG^. This model gives a point 

(r, m) = (— 71, n — 1) (45) 

rather than a curve in the m-r plane, see Fig. [l] Since r = Xi/x2 = —n, we have xi = —nx2- 
So we may consider this system in the 4 dimensional subspace X2-x^. In this subspace, while 
L2,iv3,L4 remain critical lines, Li shrinks to a critical point (1/(1 — ra), 0, 0, 1). In order to 
have the standard matter and radiation points, the point {—n,n — 1) should be located at 
(—1/2, —1/2), thus n = 1/2. Note that this result is consistent with what has been obtained 
in surveying the scaling solutions of f{G) models |10]. So in this case, all the 4 critical 
manifolds (points or lines) are at (-1/2,-1/2). Re-calculating the eigenvalues, we find 
that, for Pi, they are —4, —3, 0, —3; for P2, P3 and P4, the eigenvalues containing m'(— 1/2) 
disappear, compared to the 5 dimensional case. The eigenvalue for Pi here is not related 
to any critical line, thus the stability of this isolated non-hyperbolic point may be obtained 
by the centre manifold method [45j. Since the eigenvalue containing m'(— 1/2) for P2 has 
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In (a/a.) In (a/a.) 

3 2 

FIG. 2: The cosmic evolution for the case f{G) = a(G4 — P)s with initial conditions xi = 
-0.0025 X (1 - 10-1^), X2 = 0.005, = -0.01, X4 = 0.99951 and X5 = -0.99. Oj is the initial 
value of the scalar factor. Note that we have to set r very close to —1/2 to obtain a long enough 
matter dominated period. When oscillating into the de Sitter period, wde sts well as We// go well 
below the w = —1 divide. 



vanished, this point is either totally stable (the system will stay in the matter dominated 
era for ever and not evolve to the dark energy dominated era.) or extremely unstable (the 
system can not stay in the matter dominated era for long enough to allow cosmic structure 
to form.). Therefore this model is ruled out unless we allow extreme fine-tuning of initial 
conditions. 

The power-law model may be generalised to f{G) = a^G^ — (3y. The corresponding 
curve m(r) and its first derivative m'(r) for this model are respectively 

m(r) = -\- p — 1 and m'(r) = ^ . (46) 



To have a standard matter point, i.e., to employ the relation (31), we find it leads to a 
constraint equation for the parameters 

2pq = l. (47) 

The system needs to be able to leave the matter point, which leads to a further constraint 

g > . (48) 

Thus if these constraints are satisfied, cosmologically viable trajectories can in principle be 
obtained. 

To be specific, we set p = 3/4 and g = 2/3 (see Fig. [T| and numerically solve this model 
to get cosmic trajectories. In Fig. |2j we obtain a trajectory that is ended with a stable 
de Sitter era. For this trajectory, while the behaviours of fim, ^de and Weff are quite 
similar to those of the ACDM model in the radiation and matter dominated epochs, Weff 
oscillates rapidly and goes deep into the phantom-like region as the universe enters the de 
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FIG. 3: The cosmic evolution for the case f{G) = a(G4 — /5)3 with initial conditions xi 
-0.0025 X (1 - 5 • 10-1^), X2 = 0.005, X3 = -0.0099, X4 = 0.99951 and X5 = -0.99. 



Sitter period. On the other hand, wde at first mimics the background fluids (although 
oscillating significantly when leaving the radiation period), then at the end of the matter 
dominated epoch plunging below the w = —1 divide for a while, thereafter increasing above 
the divide, and finally oscillating rapidly into the de Sitter period. In the m-r plane, confined 
on the straight line m = 1/2 • r — 1/4, the universe slowly leaves the point (—1/2, —1/2) 
towards the r = line, bounces back and forth and then goes into the stable de Sitter point 
with r < 0. The reason why wde as well as Weff oscillate rapidly when the system enters 
the de Sitter period is related to the fact that the line r = actually corresponds to the 
bifurcation point of Li, similar to the bifurcation point of L2 that is discussed in details in 



Section III To see this, we express the last two eigenvalues of Pi in terms of r. 




(49) 



and we find that these eigenvalues have divergent imaginary parts if r ^ 0~. 

Note that the system may also evolve to the stable phantom-like point P3 with 0:50 > 1 
after the matter dominated period, see Fig. |3| In this case, wde oscillates rapidly in the 
radiation and matter dominated epochs, then plunges well below the w = —1 divide at 
the end of the matter dominated epoch and thereafter shortly passes nearby the de Sitter 
point Pi (with rapid oscillation) and the point P3 with X50 — before going into the stable 
phantom-like point. In the m-r plane's view, the system slowly leaves the point (—1/2, —1/2) 
towards the r = line and then bounces back and falls onto the point (—1/2,-1/2). 

We have also surveyed several other cases. The forms of these models, the corresponding 
m(r) (see Fig. [T]) and the constraint of the parameters from Eq. (31) are listed below: 



/(G) = aGP + pC^ 
f{G) =aGPexp{f]G) 
f{G) = aGPexp{(3/G) 
f{G) = aG^ [H(^G)Y 



m{r) = p + q — 1 + pq/r , p 

m{r) = —r + p/r , p 

— r — p/r — 2 , p 

m{r) = A{r)/qr , p 



m(r) 



1/2 kq< 1/2, 

1/2, 

1/2, 

1/2, 
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where A{r) = (1 — g)r^ + {2p — q)r + p^. The model f{G) = aG^ exp{PG) is ruled out 
(if extreme fine-tuning of initial conditions is not allowed) because the condition ( [40| can 
not be satisfied, i.e., the corresponding m(r) curve approaches the point (—1/2, —1/2) from 
a forbidden direction. For the other models, it is easy to produce an elongated matter 
dominated epoch followed by a radiation dominated epoch, but we find it difficult for the 
cosmic trajectory to evolve to a final stable accelerated epoch after the standard matter 
dominated epoch. 



V. CONCLUSIONS 

In this paper we have derived conditions for cosmologically viable f{G) models. By 
recasting the field equations of f{G) gravity with a fiat Robertson- Walker metric into a 
5 dimensional autonomous system, we have seen that the equilibria of the system are not 
isolated but form so-called equilibrium manifolds (in the present case they are 1 dimensional 
lines, called critical lines.). We have shown that the emergence of the critical lines is not 
totally due to the fact that the 5 dimensional autonomous system can in principle be reduced 
to a 4 dimensional one. Discovering 4 critical lines in the phase space for a general f{G) 
model, the curved line L3 is found to connect to the other 3 straight lines, Li, L2 and L4. 
We have also found that for all 4 lines the stability of a critical point changes when moving 
along the critical line, which may be referred to as bifurcation without parameters. For 
the critical line L2 (or L4) in which the dark energy density scales with that of matter (or 
radiation), matter (or radiation) becomes dominant only if the parameter of the critical 
line X30 is around 0. At the same time, the stability of the critical line changes rapidly 
around X30 = 0, and we require 0:30 < to avoid a violent instability (i.e., to avoid extreme 
fine-tuning of the initial conditions). 

Similar to the case of f{R), for a specific f{G) model, we can construct the m(r) curve 
in the m-r plane. In this plane, L2, and L4 all lie at the point (—1/2, —1/2). Thus for 
a standard matter point P2 (on L2) to exist, we require the m(r) curve satisfy the relation 



m — = — 
^ 2' 2 



The condition (40): 

m'i-l) > -1 

is also employed to let the universe finally leave the matter dominated point, which corre- 
sponds to some forbidden directions around the point (—1/2, —1/2) in the m-r plane. 

The stable condition for the de Sitter point Pi (on Li) corresponds to the area enclosed 
by the curve mdsir) = (r + l)/2r and the r axis in the m-r plane, in which the area enclosed 
by the curve miis*{r) = 8(r + l)/25r and the r axis corresponds to a stable spiral (in a 2 
dimensional subspace) with the rest corresponding to a stable node. P3 (on L3) is a dark 
energy dominated point. It has a wealth of possibilities. Besides mimicking radiation when 
X50 = —1 and matter when X50 = —1/2, it can be a stable quintessence-like point or a stable 



phantom-like point. But if we require a standard matter epoch (requiring the condition (40) 
to hold unless we allow extreme fine-tuning of initial conditions), it can only be stable as a 
phantom-like point. Which stable accelerated point the universe finally evolves to, the de 
Sitter point Pi or the phantom-like point P3, depends on the initial conditions. 



14 



We have also studied several toy models of f{G) dark energy whose m(r) curves can be 
obtained analytically. Note that if the analytical form of m(r) can not be obtained, one may 
obtain it numerically. On the other hand, if a form of m(r) can be obtained from observations 
or other kinds of analyses, it may not be easy to derive the corresponding f{G), but still we 
might find it numerically. For the model f{G) = a{G^ — (3Y in particular, cosmologically 
viable trajectories have been obtained for the case with a de Sitter epoch as the final stage as 
well as for the case with a phantom-like epoch as the final stage. These trajectories mimic the 
ACDM scenario in the radiation and matter dominated periods (although wde may oscillate 
for a certain period) but have distinctive signatures at late times. Most significantly, wde 
plunges below the w = —1 divide in the late matter dominated epoch (pole-like behaviour, 
for similar behaviours of wde in other scenarios see [50]) and then oscillates rapidly, also 
going below the w = —1 divide, when evolving to (for the case with a de Sitter epoch as the 
final stage) or passing nearby (for the case with a phantom-like epoch as the final stage) the 
de Sitter point Pi. For other potentially viable models, cosmologically viable trajectories 
may still exist, but we find it is non-trivial to obtain them. We leave this task for future 
work. 
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